// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                        avtCurveConstructorFilter.C                        //
// ************************************************************************* //

#include <avtCurveConstructorFilter.h>
#include <avtDataTree.h>

#include <vtkCellArray.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkVisItUtility.h>
#include <InvalidDimensionsException.h>
#include <NoInputException.h>
#include <MapNode.h>

#include <DebugStream.h>

#ifdef PARALLEL
#include <mpi.h>
#include <avtParallel.h>
#include <vtkDataSetReader.h>
#include <vtkDataSetWriter.h>
#include <vtkCharArray.h>
#endif


#include <maptypes.h>

using std::string;


// ****************************************************************************
//  Method: avtCurveConstructorFilter constructor
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Sat Apr 20 13:01:58 PST 2002
//
//  Modifications:
//    Kathleen Bonnell, Thu Mar 19 17:42:48 PDT 2009
//    Added forceConstruct.
//
//    Kathleen Bonnell, Mon Mar 23 09:53:41 PDT 2009
//    Removed forceConstruct.
//
// ****************************************************************************

avtCurveConstructorFilter::avtCurveConstructorFilter()
{
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter destructor
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Sat Apr 20 13:01:58 PST 2002
//
// ****************************************************************************

avtCurveConstructorFilter::~avtCurveConstructorFilter()
{
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter::ExecuteData
//
//  Purpose:
//      Does the actual VTK code to modify the dataset.
//
//  Arguments:
//      inDS      The input dataset.
//      <unused>  The domain number.
//
//  Returns:      The output dataset.
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Sat Apr 20 13:01:58 PST 2002
//
//  Modifications:
//
//    Hank Childs, Tue May 28 11:41:21 PDT 2002
//    Use the variable name as the label for our output data tree.
//
//    Kathleen Bonnell, Fri Jul 12 16:53:11 PDT 2002
//    Removed vtk filters associated with label-creation.  Now handled by
//    the plot.
//
//    Kathleen Bonnell, Tue Dec 23 10:18:06 PST 2003
//    Added logic to handle point-sorting when necessary. Add vertex cells,
//    so they can be displayed upon user request.
//
//    Mark C. Miller, Wed Jun  9 21:50:12 PDT 2004
//    Eliminated use of MPI_ANY_TAG and modified to use GetUniqueMessageTags
//
//    Hank Childs, Wed Jan  4 11:21:59 PST 2006
//    Allocate the size of the verts array front-end.
//
//    Kathleen Bonnell, Tue Jun 20 16:02:38 PDT 2006
//    Save the curve points in output array to be added to PlotInfoAttributes.
//
//    Kathleen Bonnell, Mon Jul 31 16:50:55 PDT 2006
//    Changed reader from PolyData to DataSetReader, as curves can now be
//    represented as 1D RectilinearGrids.  Modified logic to handle
//    Rectilinear Grids as needed.
//
//    Hank Childs, Thu Jan  4 09:29:59 PST 2007
//    Add handling for error case where a tree is not empty, but it has no
//    leaves (so it really is "IsEmpty").
//
//    Mark C. Miller, Mon Jan 22 22:09:01 PST 2007
//    Changed MPI_COMM_WORLD to VISIT_MPI_COMM
//
//    Dave Bremer, Thu Jun  7 19:47:35 PDT 2007
//    Added a fix to correctly construct a curve from a histogram plot.
//    I determine if there's a problem by looking at the geometry, and
//    if necessary flip the last two points in each group of 4.
//
//    Hank Childs, Fri Feb 15 15:52:46 PST 2008
//    Throw an exception in an error condition.
//
//    Mark C. Miller, Wed Jun 11 12:09:30 PDT 2008
//    Replaced Exception with continue
//
//    Kathleen Bonnell, Thu Mar 19 17:42:48 PDT 2009
//    Don't use unique-leaf test for early termination if forceConstruct is set.
//
//    Kathleen Bonnell, Mon Mar 23 09:54:09 PDT 2009
//    Removed unique-leaf test and forceConstruct.
//
//    Brad Whitlock, Tue Mar 30 15:33:19 PDT 2010
//    Pass along avtCurveTransform in the field data if it exists.
//
//    Kathleen Bonnell, Tue Dec 14 12:57:14 PST 2010
//    Ensure the output variable is named.
//
//    Kathleen Bonnell, Thu Feb 17 09:19:06 PST 2011
//    Moved bulk of dataset construction into CreateSingleOutput method,
//    which will allow creation of multiple outputs when necessary.
//
//    Kathleen Biagas, Thu Sep 29 06:11:38 PDT 2011
//    Only construct multiple outputs when requested in data attributes.
//
//    Kathleen Biagas, Tue Dec 19, 2023
//    Consolidate logic for no-labels and labels.
//    Add varname and count to CreateSingleOutput.

// ****************************************************************************

void avtCurveConstructorFilter::Execute()
{
    avtDataTree_p inTree = GetInputDataTree();

#ifdef PARALLEL
    //
    //  Gather all data onto one processor, so that
    //  there will be no discontinuities in the curve.
    //
    int myRank, numProcs;

    MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
    MPI_Comm_rank(VISIT_MPI_COMM, &myRank);

    int tags[3];
    GetUniqueMessageTags(tags, 3);
    int mpiNdsTag  = tags[0];
    int mpiSizeTag = tags[1];
    int mpiDataTag = tags[2];

    if (myRank == 0)
    {
        int i , j;
        for (i = 1; i < numProcs; i++)
        {
           MPI_Status stat;
           MPI_Status stat2;
           int nds = 0, size = 0;
           MPI_Recv(&nds, 1, MPI_INT, MPI_ANY_SOURCE, mpiNdsTag,
                    VISIT_MPI_COMM, &stat);
           for (j = 0; j < nds; j++)
           {
               vtkDataSetReader *reader = vtkDataSetReader::New();
               reader->ReadFromInputStringOn();
               MPI_Recv(&size, 1, MPI_INT, stat.MPI_SOURCE, mpiSizeTag,
                         VISIT_MPI_COMM, &stat2);
               char *str = new char[size];
               MPI_Recv(str, size, MPI_CHAR, stat.MPI_SOURCE, mpiDataTag,
                        VISIT_MPI_COMM, &stat2);
               vtkCharArray *charArray = vtkCharArray::New();
               charArray->SetArray((char*)str, size, 1);
               reader->SetInputArray(charArray);
               reader->Update();

               inTree->Merge(new avtDataTree(reader->GetOutput(), -1));
               delete [] str;
               reader->Delete();
               charArray->Delete();
           }
        }
    }
    else
    {
        SetOutputDataTree(new avtDataTree());
        char *str = NULL;;
        int i = 0, size = 0, nleaves = 0;
        if (inTree->IsEmpty())
        {
            MPI_Send(&nleaves, 1, MPI_INT, 0, mpiNdsTag, VISIT_MPI_COMM);
            return;
        }

        vtkDataSet **ds = inTree->GetAllLeaves(nleaves);
        MPI_Send(&nleaves, 1, MPI_INT, 0, mpiNdsTag, VISIT_MPI_COMM);

        vtkDataSetWriter *writer = vtkDataSetWriter::New();
        writer->WriteToOutputStringOn();
        writer->SetFileTypeToBinary();

        for (i = 0; i < nleaves; i++)
        {
            writer->SetInputData(ds[i]);
            writer->Write();
            size =  writer->GetOutputStringLength();
            str  =  writer->RegisterAndGetOutputString();

            MPI_Send(&size, 1, MPI_INT, 0, mpiSizeTag, VISIT_MPI_COMM);
            MPI_Send(str, size, MPI_CHAR, 0, mpiDataTag, VISIT_MPI_COMM);
            delete [] str; //allocated by writer
        }
        writer->Delete();
        delete [] ds;
        return;
    }
#endif

    if (inTree->IsEmpty())
    {
        SetOutputDataTree(inTree);
        return;
    }

    avtDataTree_p outTree;
    stringVector labels;
    inTree->GetAllLabels(labels);

    const char *vname = (pipelineVariable != NULL ? pipelineVariable : "");
    size_t count = labels.empty() ? 1 : labels.size();
    vtkDataSet **ds = new vtkDataSet *[count];
    for (size_t i = 0; i < count; ++i)
    {
        ds[i] = nullptr;
        if(labels.empty())
        {
            ds[i] = CreateSingleOutput(inTree, vname, count);
        }
        else
        {
            avtDataTree_p oneTree = inTree->PruneTree(labels[i]);
            ds[i] = CreateSingleOutput(oneTree, labels[i], count);
        }
        if (ds[i] == NULL)
        {
            for (size_t j = 0; j < i; ++j)
            {
                if (ds[j] != NULL)
                    ds[j]->Delete();
            }
            delete [] ds;
            SetOutputDataTree(inTree);
            return;
        }
    }
    if(labels.empty())
    {
       avtDataRepresentation dr(ds[0], -1, vname);
       outTree = new avtDataTree(dr);
    }
    else
    {
        outTree = new avtDataTree((int)labels.size(), ds, -1, labels);
    }

    SetOutputDataTree(outTree);
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter::CreateSingleOutput
//
//  Notes:  Moved from Execute mthod.
//
//  Purpose:
//      Does the actual VTK code to modify the dataset.
//
//  Arguments:
//      inTree    The input data tree.
//
//  Returns:      The output dataset.
//
//  Programmer: Kathleen Bonnell
//  Creation:   February 17, 2011
//
//  Modifications:
//    Kathleen Biagas, Tue Dec 19, 2023
//    Add vname and count arguments.
//    Store xy-pairs in MapNode when creating multiple curves.
//
//    Kathleen Biagas, Wed Sep 11, 2024
//    Moved test for ouputArray size from PostExecute to here, so that we
//    don't save multiple curves when outputArray is too large.
//    Utilize same varname for single-curve output as would be used for
//    multi-curve output.
//
// ****************************************************************************

vtkDataSet *
avtCurveConstructorFilter::CreateSingleOutput(avtDataTree_p inTree, const string &label, const int count)
{
    //
    //  This filter doesn't do much right now.  Basically just a
    //  "connect-the-dots" between the vertices.
    //
    int nleaves, j, k;
    double x;
    vtkDataSet **ds;
    ds = inTree->GetAllLeaves(nleaves);
    if (nleaves == 0)
    {
        // Kind of a bizarre error, since "IsEmpty" above should have returned
        // true.  This situation does occur, and it is a little dangerous to
        // change the behavior of IsEmpty, so just accomodate.
        delete [] ds;
        return NULL;
    }
    DoubleIntMap minX;

    //
    //  Cannot assume that there are no overlaps of points
    //  between datasets, so keep track of first and last
    //  x values, to determine if sorting is required.
    //  Now must order the datasets so that the points will
    //  be accessed in proper order, to avoid discontinuities
    //  between domains.
    //  Store in map with x value as the key.
    //  (from first point in each ds).
    //
    double *mm = new double[nleaves*2];
    int npts;

    for (j = 0; j < nleaves; j++)
    {
        vtkDataArray *data = NULL;
        if (ds[j]->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        {
            data = ((vtkRectilinearGrid*)ds[j])->GetXCoordinates();
        }
        else if (ds[j]->GetDataObjectType() == VTK_POLY_DATA)
        {
            data = ((vtkPolyData*)ds[j])->GetPoints()->GetData();
        }
        else
        {
            continue;
        }
        npts = data->GetNumberOfTuples();
        x = data->GetComponent(0, 0);
        minX.insert(DoubleIntMap::value_type(x, j));
        mm[j*2] = x;
        mm[j*2+1] = data->GetComponent(npts-1, 0);
    }

    bool requiresSort = false;
    for (j = 0; j < nleaves && !requiresSort; j++)
    {
        for (k = 0; k < nleaves && !requiresSort; k++)
        {
            if (k != j && mm[k*2] > mm[j*2] && mm[k*2] < mm[j*2+1])
                requiresSort = true;
        }
    }
    delete [] mm;

    vtkDataArray *inXC;
    vtkDataArray *inVal;

    int dtype = VTK_FLOAT;
    if (ds[0]->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        dtype = ((vtkRectilinearGrid*)ds[0])->GetXCoordinates()->GetDataType();
    else
        dtype = ((vtkPolyData*)ds[0])->GetPoints()->GetData()->GetDataType();
    vtkRectilinearGrid *outGrid = vtkVisItUtility::Create1DRGrid(0, dtype);
    vtkDataArray *outXC  = outGrid->GetXCoordinates();
    vtkDataArray *outVal = outXC->NewInstance();
    outGrid->GetPointData()->SetScalars(outVal);
    outVal->Delete();
    int nPoints;

    //
    //  Ensure that the output is 2d by setting z-component to zero.
    //
    DoubleIntMap::iterator it;
    int index = 0;
    for (it = minX.begin(); it != minX.end(); it++)
    {
        if (ds[(*it).second]->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        {
            inXC  = ((vtkRectilinearGrid*)ds[(*it).second])->GetXCoordinates();
            inVal = ((vtkRectilinearGrid*)ds[(*it).second])->GetPointData()->
                    GetScalars();
            nPoints = inXC->GetNumberOfTuples();
            for (vtkIdType i = 0; i < nPoints; i++, index++)
            {
                outXC->InsertNextTuple1(inXC->GetTuple1(i));
                outVal->InsertNextTuple1(inVal->GetTuple1(i));
            }
        }
        else if (ds[(*it).second]->GetDataObjectType() == VTK_POLY_DATA)
        {
            inXC  = ((vtkPolyData*)ds[(*it).second])->GetPoints()->GetData();
            nPoints = inXC->GetNumberOfTuples();

            // Insert a check to see if this data comes from a histogram plot,
            // in which case we need to swap some of the points.  Data comes
            // in groups of 4 points for each bar, connected like this: N
            // so we reverse the last two points to create bar shapes instead.
            bool bIsHistogramPlot = false;
            if (nPoints%4 == 0)
            {
                bIsHistogramPlot = true;
                for (int i = 0; i < nPoints; i+=4)
                {
                    if ( inXC->GetComponent(i,   0) != inXC->GetComponent(i+1, 0) ||
                         inXC->GetComponent(i+2, 0) != inXC->GetComponent(i+3, 0) ||
                         inXC->GetComponent(i,   1) != inXC->GetComponent(i+2, 1) ||
                         inXC->GetComponent(i+1, 1) != inXC->GetComponent(i+3, 1) )
                    {
                        bIsHistogramPlot = false;
                        break;
                    }
                }
            }
            if (!bIsHistogramPlot)
            {
                for (int i = 0; i < nPoints; i++, index++)
                {
                    outXC->InsertNextTuple1(inXC->GetComponent(i, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i, 1));
                }
            }
            else
            {
                for (int i = 0; i < nPoints; i+=4, index+=4)
                {
                    outXC->InsertNextTuple1( inXC->GetComponent(i,   0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i,   1));
                    outXC->InsertNextTuple1( inXC->GetComponent(i+1, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i+1, 1));
                    outXC->InsertNextTuple1( inXC->GetComponent(i+3, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i+3, 1));
                    outXC->InsertNextTuple1( inXC->GetComponent(i+2, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i+2, 1));
                }
            }
        }
    }

    //
    // Sort if necessary
    //
    vtkDataArray *sortedXC;
    vtkDataArray *sortedVal;
    if (requiresSort)
    {
        sortedXC = outXC->NewInstance();
        sortedVal = outVal->NewInstance();
        DoubleIntMap sortedIds;
        nPoints = outXC->GetNumberOfTuples();
        for (int i = 0; i < nPoints; i++)
        {
            x = outXC->GetTuple1(i);
            sortedIds.insert(DoubleIntMap::value_type(x, i));
        }
        DoubleIntMap::iterator it;
        for (it = sortedIds.begin(); it != sortedIds.end(); it++)
        {
            sortedXC->InsertNextTuple1(outXC->GetTuple1((*it).second));
            sortedVal->InsertNextTuple1(outVal->GetTuple1((*it).second));
        }
        outGrid->SetXCoordinates(sortedXC);
        outGrid->GetPointData()->SetScalars(sortedVal);
        sortedXC->Delete();
        sortedVal->Delete();
    }
    else
    {
        sortedXC  = outXC;
        sortedVal = outVal;
    }
    nPoints = sortedXC->GetNumberOfTuples();
    outGrid->SetDimensions(nPoints, 1, 1);
    outputArray.clear();
    for (int i = 0; i < nPoints; i++)
    {
        outputArray.push_back(sortedXC->GetTuple1(i));
        outputArray.push_back(sortedVal->GetTuple1(i));
    }
    string varname = (!label.empty() ? label : (pipelineVariable != NULL ? pipelineVariable : "Curve"));
    // Limit outputArray size that we send back to the client.
    if(outputArray.size() < 100000)
    {
        if (count > 1)
        {
            outputInfo[varname] = outputArray;
            outputArray.clear();
        }
    }
    else
    {
        debug5 << "Curve constructor filter does not send curves that contain "
                  "more than 100K values to the client." << endl;
    }

    // Pass a copy of the avtCurveTransform data into the new curve object.
    vtkDataArray *ct = ds[0]->GetFieldData()->GetArray("avtCurveTransform");
    if(ct != 0 && ct->GetNumberOfTuples() == 16)
    {
        ct->Register(NULL);
        outGrid->GetFieldData()->AddArray(ct);
    }


    // make sure the outputvar is named.
    sortedVal->SetName(varname.c_str());


    //
    //  Clean up.
    //
    delete [] ds;
    return outGrid;
}

// ****************************************************************************
//  Method: avtCurveConstructorFilter::VerifyInput
//
//  Purpose:
//      Verifies that the input is 2D data, throws an exception if not.
//
//  Programmer: Kathleen Bonnell
//  Creation:   April 26, 2002
//
// ****************************************************************************

void
avtCurveConstructorFilter::VerifyInput(void)
{
    if  (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() != 1)
    {
        EXCEPTION2(InvalidDimensionsException, "Curve", " Lines ");
    }
}

// ****************************************************************************
//  Method: avtCurveConstructorFilter::ModifyContract
//
//  Purpose:
//    Indicates that we cannot do dynamic load balancing with this filter.
//
//  Programmer: Kathleen Bonnell
//  Creation:   April 26, 2002
//
//  Modifications:
//
//    Hank Childs, Tue Feb 19 19:45:43 PST 2008
//    Rename "dynamic" to "streaming", since we really care about whether we
//    are streaming, not about whether we are doing dynamic load balancing.
//    And the two are no longer synonymous.
//
// ****************************************************************************

avtContract_p
avtCurveConstructorFilter::ModifyContract(avtContract_p spec)
{
    spec->NoStreaming();
    return spec;
}

// ****************************************************************************
//  Method: avtCurveConstructorFilter::ModifyContract
//
//  Purpose:
//      Allows the filter to change its output's data object information, which
//      is a description of the data object.
//
//  Programmer: Kathleen Bonnell
//  Creation:   December 23, 2002
//
//  Modifications:
//    Kathleen Bonnell, Mon Jul 31 16:50:55 PDT 2006
//    Changed Spatial dimension to 1, as curves are not represented as
//    1D RectilinearGrid.
//
// ****************************************************************************

void
avtCurveConstructorFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetAttributes().SetSpatialDimension(1);
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter::PostExecute
//
//  Purpose:
//    Saves the curve data to an output array in PlotInfoAtts.
//
//  Programmer: Kathleen Bonnell
//  Creation:   June 20, 2006
//
//  Modifications:
//    Brad Whitlock, Tue Jan  6 16:41:18 PST 2009
//    Store the outputArray into the MapNode data.
//
//    Brad Whitlock, Mon Jul 23 12:15:34 PDT 2012
//    Limit the number of values sent to the client in plot information.
//
//    Kathleen Biagas, Tue Dec 19, 2023
//    Multi-curve output is stored in outputInfo. Single curve stored in
//    outputArray.  Removed parallel broadcast (don't think its necessary).
//
//    Kathleen Biagas, Wed Sep 11, 2024
//    Moved test for ouputArray size to CreateSingleOutput, so that we
//    don't save multiple curves when outputArray is too large.
//
// ****************************************************************************

void
avtCurveConstructorFilter::PostExecute(void)
{
    if(!outputArray.empty())
    {
        PlotInfoAttributes plotInfoAtts;
        MapNode data;
        data = outputArray;
        GetOutput()->GetInfo().GetAttributes().AddPlotInformation("Curve", data);
    }
    else if (outputInfo.GetNumEntries() > 0)
    {
        // There were multiple curves
        GetOutput()->GetInfo().GetAttributes().AddPlotInformation("Curves", outputInfo);
    }
}

